### Expected Utility Maximization with Risk Aversion function library
## Appendix Theory: A.4.1

MaxLikeSR_CRRA_LOO_crossVal = function(data, sub_sample, init){
  # Leave-one-out Cross Validation (LOO CV)
  start <- Sys.time()  
  data = sample_SR(data,sub_sample)
  index = data$studentID                                                            
  data = split(data,f = index)  
  
  dfs_LOOCV = lapply(data, LOOCV_train_test_split)
  results = lapply(dfs_LOOCV, LOO_crossVal_train_predict_CRRA, init)
  CRRA_coef = unlist(lapply(results, {function(x) mean(x[["CRRA.coef"]])}), use.names = FALSE)
  CRRA_coef_std = unlist(lapply(results, {function(x) sqrt(var(x[["CRRA.coef"]]))}), use.names = FALSE)
  score = unlist(lapply(results, {function(x) sum(x$accuracy)}), use.names = FALSE)
  results_scores = as.data.frame(cbind(studentID = unique(index), score = score, CRRA_mean = CRRA_coef,
                                       CRRA_std = CRRA_coef_std))
  
  print("Constant Relative Risk Aversion - Leave-One-Out Test Accuracy Summary:")
  print(paste("Average:", mean(results_scores$score)))
  print(paste("Standard error:", sqrt(var(results_scores$score)/length(results_scores$score))))
  
  end = Sys.time()                                                                  
  timer = difftime(end,start,units = "secs")
  print(paste("Time taken for MLE algorithim and Cross-Validation:",seconds_to_period(round(timer[[1]],2))))
  return(list(estimations = results, summary = results_scores))
}



LOOCV_train_test_split = function(data){
  result = list()
  for(i in 1:length(data$qnum)){
    test_df = data[i, ]
    train_df = data[-i, ]
    result[[i]] <- list(train = train_df, test = test_df)
  }
  result
}



LOO_crossVal_train_predict_CRRA = function(LOO_train_test_folds, init){
  results = bind_rows(lapply(LOO_train_test_folds, LOO_train_test_CRRA, init))
  results
}



LOO_train_test_CRRA = function(train_test_dfs, init){
  train_df = train_test_dfs[["train"]]
  test_df = train_test_dfs[["test"]]
  
  results = wrapMLE_CRRA(train_df, init, log_lik_CRRA)
  results_summary = stat_sum_CRRA(results)
  results_summary$studentID = test_df$studentID
  
  param = list(CRRA = results_summary$CRRA.coef)
  test_predict_results = LOO_test_predict(test_df, param, calc_util_CRRA)
  results_summary = cbind(results_summary, test_qnum = test_df$qnum, test_choice = test_df$choice, test_predict_results)
  results_summary
}





LOO_test_predict = function(test_df, param, util_fun){
  lottery_data = list(l1 = list(prob = test_df[["prob.1"]], lower = test_df[["lower.1"]], upper = test_df[["upper.1"]]),
                      l2 = list(prob = test_df[["prob.2"]], lower = test_df[["lower.2"]], upper = test_df[["upper.2"]]),
                      l3 = list(prob = test_df[["prob.3"]], lower = test_df[["lower.3"]], upper = test_df[["upper.3"]]),
                      l4 = list(prob = test_df[["prob.4"]], lower = test_df[["lower.4"]], upper = test_df[["upper.4"]]),
                      l5 = list(prob = test_df[["prob.5"]], lower = test_df[["lower.5"]], upper = test_df[["upper.5"]]))
  
  lottery_utility = unlist(lapply(lottery_data, util_fun, param))
  predicted_choice = unname(which.max(lottery_utility))
  accuracy = if_else(predicted_choice == test_df$choice, 1, 0)
  list(pred_choice = predicted_choice, accuracy = accuracy)
}


calc_util_CRRA = function(df, param){
  if(param[["CRRA"]]<=1.01 & param[["CRRA"]]>=0.99){
    df[["prob"]]*log(55+df[["upper"]]) + (1-df[["prob"]])*log(55+df[["lower"]])
  }
  else{
    df[["prob"]]*(((55+df[["upper"]])^(1-param[["CRRA"]]))/(1-param[["CRRA"]])) + (1-df[["prob"]])*(((55+df[["lower"]])^(1-param[["CRRA"]]))/(1-param[["CRRA"]]))
  }
}



wrapMLE_CRRA = function(data,init,logLikFun){
  CRRA_MaxLike_res = mle2(minuslogl = logLikFun, data = list(data=data), start = init, 
                         lower=c(CRRA=0))
  CRRA_MaxLike_res
}


stat_sum_CRRA = function(stat){
  results = data.frame("studentID"=NA,"CRRA.coef"=NA,
                       "std.CRRA"=NA,"p.val.CRRA"=NA,
                       "logLikelihood"=NA,"AIC"=NA)
  summary = summary(stat)
  results["studentID"]  = NA
  results["CRRA.coef"]  = summary@coef[1,1]
  results["std.CRRA"]   = summary@coef[1,2]
  results["p.val.CRRA"] = summary@coef[1,4]
  results["logLikelihood"] = logLik(stat)
  results["AIC"]  = AIC(stat)
  return(results)
}



log_lik_CRRA <- function(CRRA, data){
  if(CRRA<=1.01 & CRRA>=0.99){
    #calculating denominator of likelihood of a choice
    denomCalc = function(data){
      exp(data[["prob.1"]]*log(55+data[["upper.1"]]) + (1-data[["prob.1"]])*log(55+data[["lower.1"]])) +
      exp(data[["prob.2"]]*log(55+data[["upper.2"]]) + (1-data[["prob.2"]])*log(55+data[["lower.2"]])) +
      exp(data[["prob.3"]]*log(55+data[["upper.3"]]) + (1-data[["prob.3"]])*log(55+data[["lower.3"]])) +
      exp(data[["prob.4"]]*log(55+data[["upper.4"]]) + (1-data[["prob.4"]])*log(55+data[["lower.4"]])) +
      exp(data[["prob.5"]]*log(55+data[["upper.5"]]) + (1-data[["prob.5"]])*log(55+data[["lower.5"]]))
    }
    
    denomRes = apply(data, 1, denomCalc) #vector of denominator components for each student
    
    # calculate numerator of likelihood of a choice
    numerCalc = function(data){
      j = data[["choice"]]
      numercalc = exp(data[[13+j]]*log(55+data[[8+j]]) + (1-data[[13+j]])*log(55+data[[3+j]]))
      numercalc
    }
    
    numerRes = apply(data, 1, numerCalc)
    
    #calculating log likelihood
    log_lik_i <-  log(numerRes/denomRes)
    
    -sum(unlist(log_lik_i))
    
  }
  else{
    denomCalc = function(data){
      exp(data[["prob.1"]]*(((55+data[["upper.1"]])^(1-CRRA))/(1-CRRA)) + (1-data[["prob.1"]])*(((55+data[["lower.1"]])^(1-CRRA))/(1-CRRA))) +
      exp(data[["prob.2"]]*(((55+data[["upper.2"]])^(1-CRRA))/(1-CRRA)) + (1-data[["prob.2"]])*(((55+data[["lower.2"]])^(1-CRRA))/(1-CRRA))) +
      exp(data[["prob.3"]]*(((55+data[["upper.3"]])^(1-CRRA))/(1-CRRA)) + (1-data[["prob.3"]])*(((55+data[["lower.3"]])^(1-CRRA))/(1-CRRA))) +
      exp(data[["prob.4"]]*(((55+data[["upper.4"]])^(1-CRRA))/(1-CRRA)) + (1-data[["prob.4"]])*(((55+data[["lower.4"]])^(1-CRRA))/(1-CRRA))) +
      exp(data[["prob.5"]]*(((55+data[["upper.5"]])^(1-CRRA))/(1-CRRA)) + (1-data[["prob.5"]])*(((55+data[["lower.5"]])^(1-CRRA))/(1-CRRA)))
    }
    
    denomRes = apply(data, 1, denomCalc) #vector of denominator components for each student
    
    # calculate numerator of likelihood of a choice
    numerCalc = function(data){
      j = data[["choice"]]
      numercalc = exp(data[[13+j]]*(((55+data[[8+j]])^(1-CRRA))/(1-CRRA)) + (1-data[[13+j]])*(((55+data[[3+j]])^(1-CRRA))/(1-CRRA)))
      numercalc
    }
    
    numerRes = apply(data, 1, numerCalc)
    
    #calculating log likelihood
    log_lik_i <-  log(numerRes/denomRes)
    
    -sum(unlist(log_lik_i))
  }
}


sample_SR = function(data,sub_sample){
  if(sub_sample == "full"){data = data}
  if(sub_sample == "part1"){data = data %>% filter(qnum<19)}
  if(sub_sample == "part2"){data = data %>%  filter(qnum>18)}
  return(data)
}






